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Metapopulation dynamics in a complex ecological landscape 
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We propose a general model to study the interplay between spatial dispersal and environment spa- 
tiotemporal fluctuations in metapopulation dynamics. An ecological landscape of favorable patches 
is generated like a Levy dust , which allows to build a range of patterns, from dispersed to clustered 
ones. Locally, the dynamics is driven by a canonical model for the evolution of the population 
density, consisting of a logistic expression plus multiplicative noises. Spatial coupling is introduced 
by means of two spreading mechanisms: diffusive dispersion and selective migration driven by patch 
suitability. We focus on the long-time population size as a function of habitat configurations, envi¬ 
ronment fluctuations and coupling schemes. We obtain the conditions, that the spatial distribution 
of favorable patches and the coupling mechanisms must fulfill, to grant population survival. The 
fundamental phenomenon that we observe is the positive feedback between environment fluctuations 
and spatial spread preventing extinction. 


PACS numbers: 87.23.Cc, 89.75.Fb, 05.40.-a 

I. INTRODUCTION 

Habitat fragmentation is commonly observed in na¬ 
ture associated with heterogeneity in the distribution of 
resources, e.g., water, food, shelter sites, physical fac¬ 
tors such as light, temperature, moisture, and any fea¬ 
ture able to affect the growth rate of the population of a 
given species [I]. A fragmented population made of sub¬ 
populations receives in the literature the suitable name 
of metapopulation cm]. These fragments, also known as 
patches, are not completely isolated as they are coupled 
due to movements of individuals in space. For model¬ 
ing purposes, as a first step one can adopt a single patch 
viewpoint, taking into account the impact of the sur¬ 
rounding population in an effective manner HO- As a 
further step beyond the single patch level, one can re¬ 
sort to a spatially explicit model. From this perspec¬ 
tive, deterministic and stochastic theoretical models have 
been developed to obtain the macroscopic behavior of the 
whole population mwm- One of the main results is 
the detection of critical thresholds that delimit the con¬ 
ditions for the sustainability of the population, which oc¬ 
curs for a suitable combination of diverse factors, related 
to quality and spatial structure of the habitat, migration 
strategies and extinction rates. Here, we address related 
fundamental questions in metapopulation theory propos¬ 
ing a model that includes a general dispersion process, 
incorporating random and selective dispersal strategies. 
Additionally, we investigate the model dynamics on top 
of a complex ecological landscape whose spatial structure 
can be tuned, ranging from spread to aggregated patches. 

Let us start from the local dynamics perspective. Lo¬ 
cally, each patch has its dynamics driven primarily by 
reproduction and intraspecific competition (carrying ca- 
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pacity). Therefore, we assume that the deterministic fac¬ 
tors that rule the evolution of the local population can 
be modeled by the logistic or Verhulst expression m 
Stochasticity is introduced in real systems by the inher¬ 
ent complexity of the fluctuating environment (external 
noise) or by the variations in the birth-death process (in¬ 
ternal, demographic noise ) n m m m, hence it has 
to be also taken into account. These deterministic and 
stochastic rules, along the lines of the canonical mod¬ 
eling constitute the local component of our 

model. We assume that this local dynamics takes place 
on each site of a lattice, where we construct a complex ar¬ 
rangement of favorable and unfavorable patches. We de¬ 
fine as favorable patches those sites that induce positive 
growth at low densities and as unfavorable patches those 
that are adverse to support life. A typical configuration 
of the model system in a square lattice is illustrated in 

Fig- 0 
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Figure 1. Ecological landscape (a green cross denotes the pres¬ 
ence of a favorable patch) and population density distribution 
(in gray scale). 







2 


Spatial coupling is introduced by migrations from one 
patch to another. First, let us assume that spatial spread 
is conservative, preserving the number of individuals dur¬ 
ing travels and also that it is nonlocal, in the sense that 
individuals can travel long distances over the landscape, 
for example like butterflies and birds mm- We model 
the populational exchange between patches based on two 
behavioral strategies: one where the individuals spread 
in space diffusively, driven by density differences, and an¬ 
other where individuals transit selectively, mainly driven 
by patch-quality differences. These strategies can be 
linked to the amount of spatial information acquired by 
the individuals HD. If they do not have any informa¬ 
tion about the spatial distribution of favorable patches, 
random movements emerge. In fact, this has been the fo¬ 
cus of works on animal foraging, where optimal efficiency 
in resource search occurs without previous knowledge of 
food distribution Il8l . This type of behavior has isotropy 
as a main trait, indicating directional indifference. In¬ 
stead, if individuals have information on the distribution 
of favorable patches, some directions will be preferred. 

In Sec. [TTJ we will describe in detail each part of the 
model. The spreading process and spatial configuration 
of the ecological landscape are described in Secs. |II A 
and II B[ respectively. The results, reported in Secs. |III 
and IV[ focus on the impact of the spatial arrangement 
of the habitat on its overall viability, that is, on the long¬ 
time behavior of the population size. Mainly numerically, 
and with the aid of analytical considerations, we inves¬ 
tigate the impact of habitat topology, spread range and 
stochasticity in the long time behavior of the population 
size, compared to the corresponding uncoupled metapop¬ 
ulation. 


II. MODEL 


In mathematical terms, we assume that the evolution 
of the population density (number of individuals per unit 
area), iq, in each patch i is described by 

Ui = CLiUi - bvif + DTi[u] + a v m • r/i(t ) + a^^/ul o &(i), 

( 2 . 1 ) 

where the constants a*, b, D are the growth rate, the 
intraspecific competition coefficient and the spatial cou¬ 
pling coefficient. The noises rji and introduce environ¬ 
ment and demographic stochasticities, respectively. They 
are assumed to be mutually independent zero mean and 
unit variance Gaussian white noises. The intensities of 
these multiplicative noises are controlled by coefficients 
<j v and 07 . The environmental noise term is expected to 
have external origins, then, its correlation even if small 
is non-null, justifying the use of Stratonovich calculus 
(•); at the same time, the demographic noise represents 
fluctuations in the reproduction process of each inde¬ 
pendent individual, then Ito calculus is more suitable 
(o) HHM. These contributions define a local dynam¬ 
ics, at each site, ruled by Eq. (2.1) with D = 0, which is 


known as canonical model [ T| . W e incorporate the addi¬ 
tional term DTi[u\ into Eq. (2.1l to account for the non¬ 
local contribution arising from fluxes between patches, as 
it will be explained below. This is the term that couples 
the set of stochastic differential equations ( 2 . 1 ). 

The intrinsic growth rate of each patch i is quantified 
by the growth rate cq, that can take positive or negative 
values. Patches can be favorable (or not), promoting 
growth (or decrease) of the local population, with eq = 
A+ >0 (or cq = A i < 0). For the sake of simplicity, we 
consider a binary landscape, where sites can be in any of 
two states, A+ = —A~ = A > 0, as assumed in previous 
studies [2D, !2T]. 


A. Nonlocal coupling 


In order to define the coupling scheme let us state some 
considerations. First, note that it is reasonable to assume 
that active individuals like butterflies, birds, terrestrial 
animals use their perception and memory to increase the 
efficiency in the search for viable habitats. The spatial 
information stored by the individuals can yield optimized 
routes between favorable regions. In fact, there is a rela¬ 
tion between spatial memory and migration strategy m- 
We introduce this trait by allowing individuals to have 
access to information about the spatial distribution of 
patch quality. Spatial knowledge can be acquired, for in¬ 
stance, by a direct verification in a previous visit or by the 
perception of the collective dynamics. Otherwise, if indi¬ 
viduals do not have any information about the ecological 
landscape, or if they do not have memory, uncorrelated 
trajectories can emerge. 

We contemplate both scenarios by modeling spread 
through a diffusive component together with a contribu¬ 
tion of direct routes connecting favorable patches, gov¬ 
erned by quality differences. The relative contribution 
of both mechanisms is regulated by parameter <5, with 
0 < 5 < 1 tuning from the ecologically driven (<5 = 0) to 
the purely diffusive (<5 = 1) cases. Moreover, we assume 
that coupling is weighted by a factor 7 (dij) that decays 
with the distance d^ between patches i and j, as will be 
defined below. Then, the flux from patch i to j is 
given by 

Jij = [5 + (1 - 5)aij] y( dij)ui > 0, (2.2) 

where = (aj — cq)/(4A) + 1/2. Hence, the total flux 
is 

Ti[n] = ^ ^ (Jji ~~ Jij) 
j^i 

= ^7 (dij) [5{uj - m) + (1 - S){ajiUj - a^m )] . 

(2.3) 


The tot al de nsity is conserved by the exchanges described 
by Eq. (2.3), as can be seen by summing over i. It indi¬ 
cates that individuals tend to move towards patches with 
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fewer individuals and better quality. For <5=1, Eq. (2.3) 
represents a generalization of the Fick’s law for nonlocal 
dispersal driven by density gradients. For <5 = 0, with 
our definition of a , t] , and binary patch growth rate, the 
possible values of oijjUj — a-,; ; yit,; are 


\ ^ 
j\ 

A 

-A 

A 

(■Uj -Ui)/2 

Uj 

-A 

- Ui 

(Uj - Uj)/2 


This means that, when the quality of two patches is dif¬ 
ferent, the flux occurs in the direction of the higher qual¬ 
ity, weighted by the out-flowing population density (low¬ 
est quality patch). Only when the quality is the same, 
diffusive exchange can occur, to allow a network of favor¬ 
able patches. 

Concerning the factor that takes into account the dis¬ 
tance between patches, there is empirical evidence {3 E2] 
that the frequency of occurrence of flights between 
patches decays with the distance, which is reasonable due 
to the increase of energetic cost. Although diverse decay 
laws are possible, we will assume exponential decay of 
the weight 7 with the traveled distance £, as observed for 
some kinds of butterflies [SI [221 (23 , that is 


l{£) = Af 1 exp(-£/4), (2.4) 


where £ c is a characteristic length (the average traveled 
distance) and the normalization constant Af is such that 
the sum of the contributions of all patches equals one. 
Operationally, we will truncate the exponential at t ~ 
8 £ c << L, where L is the linear characteristic size of the 
landscape. 


widely spread (for p = 0 ) to compactly aggregated in a 
few clusters separated by large empty spaces (for p = 3), 
as illustrated in Fig. [2] 
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Figure 2. Habitat topology for different values of the exponent 
p. Black cells indicate positive growth rate A and white cells 
negative growth —A, in a square domain of linear size L = 
100. The density of favorable patches is p = 0.1. 

We quantify the change in the spatial structure by 
computing the probability distribution of the minimal 


B. Ecological landscape 

In nature, the arrangement of the ecological landscape 
is built by many distinct processes, occurring in many 
time scales, creating complex spatiotemporal structures. 
Then, beyond the inclusion of the environmental noise 
77, it is also important to take into account the spatial 
organization of patches [HI EH [24J [25]. 

Heterogeneity and patchiness are adequate to cap¬ 
ture the complexity of diverse ecological systems ES¬ 
SO]. Here we propose to use as complex ecological land¬ 
scape a Levy dust [lS] distribution of favorable patches 
on a square domain of size L x L patches, with peri¬ 
odic boundary conditions. Over a background of adverse 
patches (a* = — A ), we construct a Levy dust of favor¬ 
able patches (at = A) given by the sites visited by a Levy 
random walk with step lengths £ drawn from the prob¬ 
ability density function p(£) a 1/£ M , with 1 < i < L. 
This protocol has been used in the study of different 
problems [18| [22 El] , but we apply it here in the study 
of metapopulation dynamics. It allows to mimic a gen¬ 
eral class of realistic conditions [23 [2BH5D] and to tune 
different habitat landscapes through parameter p, from 



Figure 3. Probability distribution of the minimal distance 
between favorable patches for different values of p, for p = 0.1 
and L = 100 (100 configuration were used). Fluctuations are 
due to the discrete nature of the possible distances in the 
lattice. The dotted lines are a guide to the eye. The solid 
line represents the probability distribution for the distance 
between uniformly distributed random points in continuous 
space, drawn for comparison. 
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distance d between favorable patches P M (d) (see Fig. [3]). 
For the density p = 0.1 used in the figure, when p ^1, 
patches are typically far from each other. For high val¬ 
ues (p > 3), the generating walk approaches the standard 
random walk, creating a much more clustered structure, 
evidenced by the peak at short distances. However the 
shape of P^{d) changes with p. When the patch density 
p is high, the shape of P M (d) resembles that of the uni¬ 
form arrangement even for large p , while at low densities 
Pfj,{d) presents a peak at small d since the resulting con¬ 
figuration of patches is very localized even for small p, 
as will be discussed in Sec. IV C Furthermore, P M (d) is 


also sensitive to L, but we kept L fixed (L = 100), even 
if some properties may have not attained the large size 
limit, as far as p and p allow to scan many qualitatively 
different possibilities of landscape structure. 

Concerning the fac tor that reflects the topology, 
as defined in Eq. (A3), it is affected by p more through 
the amount of favorable patches n v than by its indirect 
consequences on the spatial distribution P M . 


will conserve the initial population size, so no interest¬ 
ing phenomena occur from the viewpoint of extinction. 
However, when 5 > 0, the diffusive behavior induces ex¬ 
ploration of the neighborhood independently of habitat 
quality, which leads to the occupation of unfavorable re¬ 
gions making likely the death of individuals. 


By numerical integration of Eq. (2.1) we obtain the 


time evolution of the total population density n(t) = 

Ef=i «*(*)• In Fig. 0 we show the outcomes for fixed 
values of the model parameters and different initial con¬ 
ditions (different landscapes). While some of the realiza¬ 
tions lead to exponential decay of the population other 
ones attain finite values at long times. Several different 
non null steady states can be attained. Notice however, 
that the steady values of different realizations are all be¬ 
low that of the uncoupled case, pi? A/b — 1000 for the 
parameters of the figure. Hence, diffusion favors the de¬ 
crease of the total population density and the occurrence 
of extinctions, as expected. 


C. General considerations about the model 


The set of parameters {D,5,1c} regulate the nonlocal 
dynamics. While D is the strength of the nonlocal cou¬ 
pling, 5 controls the balance between diffusion and di¬ 
rected migration, and £ c defines the coupling range. The 
ecological landscape is characterized by p and p. 

In the results presented in the following sections, we 
will restrict the analysis to a region of parameter space 
relevant to discuss the main phenomenology of the model. 
Thus, we will set A = b = 1 in all cases. We will also 
consider L = 100 and typically p = 0.1. Concerning 
the noise parameters, we set a v = = 0 to analyze the 

deterministic case in Sec. |III| and turn noise on by setting 
<j v = <j£ = 1 in Sec. IV This choice is based on previous 
works mm- Indeed, population size can be subject 
to large fluctuations as demonstrated by experimental 
data [TUI . 


We performed numerical simulations of Eq. (2.1) on 


top of different landscapes, by preparing the system in 
the stationary state of the deterministic and uncoupled 
case, i.e., rq(0) = max{ qj/b , 0} for all i, plus a small 
noise. Integration of Eq. (2.1) was carried out with Euler- 


Maruyama scheme with a time step At = 10 
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III. DETERMINISTIC CASE (<7„ = = 0) 

Before proceeding to study the full model, we consider 
the deterministic case. Locally, when stochastic contri¬ 
butions are neglected, the asymptotic value of the popu¬ 
lation size for each patch is Ui = ai/b. Introducing non¬ 
local effects, the population size might change. If pop¬ 
ulation exchanges between patches are guided solely by 
their quality (<5 = 0), then, the favorable-patch network 



t 

Figure 4. Deterministic (a v = = 0) time evolution of 

the total population density n, for <5 = 1, p = 0.1, D = 10, 
4 = 0.5, p = 1.7, and different initial landscapes. This set 
of values results in about half of 50 realizations leading to 
extinction. We use a dotted line to flag the ones that tend 
to extinction exponentially fast and a solid line for those that 
lead to population survival. 


In order to investigate how the fraction of survivals 
changes with the topology, we plot in Fig. [5] the num¬ 
ber of survivals per realization, f s , as function of p, for 
several values of D. Besides the initial condition used 
throughout this paper (see Sec. IIC), we observed that 
a perturbation of the null state also leads to the same 
results of Fig. [5] For given p, increasing D favors the 
occurrence of extinctions as already commented above. 
For given D, below a threshold value of p the population 
gets extincted in all the realizations, while above a sec¬ 
ond threshold it always survives (for the finite number of 
realizations done), between thresholds both states, the 
null and non null ones, are accessible. The number of 
non null stable states increases with p. 
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Summing over all i the deterministic form of Eq. (2.1), 
one finds that the steady solution h = 0 must satisfy 
Aytq = bJ2 u ii which has infinite solutions between 
the fundamental null state and the uncoupled case solu¬ 
tion (the only stable one for D = 0). The condition for 
stationarity of the total density depends only on the local 
parameters, since fluxes are only internal, however, the 
coupling and landscape can stabilize configurations other 
than the trivial ones. Furthermore, in the Appendix, we 
performed an approximate calculation to show that, for 
small D , the null state is stable if 


A-D( l- 7 / 1 )>0, (3.1) 

where 0 < 7 ^ < 1 is a factor that mirrors the topol¬ 
ogy, varying from 7 ^ = p for the uniform case p = 0 
to 7 ^ = 1 in the limits of large p or large p. Despite 
this approximate expression fails in providing accurate 
threshold values, it predicts that survival is facilitated 
by larger A and spoiled by increasing D. It also qual¬ 
itatively predicts the impact of the topology, as far as 
it indicates that the destructive role of diffusion can be 
compensated by a large enough degree of clusterization 
of the resources given by large 7 M . 


IV. STOCHASTIC CASE 

First let us review some known results about the local 
(one site) dynamics, which is obtained in the limit D —» 0 
of Eq. |2.1| (canonical model). In the deterministic case, 
the two-state habitat [ 20 ] [ 2 T leads to local extinction 
(if cii = —A ) or finite population (if a* = +A). The 
presence of stochastic contributions changes the stability 
of the patches. When a* = —A < 0, the local extinction 
event predicted deterministically is reinforced by noise. 



Figure 5. Fraction of surviving metapopulations f s (over 100 
realizations) in the deterministic case (oy = ay = 0) as a 
function of p, setting 5 = 1, p = 0.1, i c = 0.5, for the values 
of D indicated on the figure. In this and following figures, 
dotted lines are a guide to the eye. 



Figure 6. Single patch dynamics. Mean extinction T time vs 
ay (for oy = 1) and vs ay (for a£ = 1). Symbols correspond 
to numerical simulations averaged over 500 sampl es an d the 
full lines to the theoretical prediction given by Eq. (4.11. The 
curve for variable ay diverges in the limit ay —> 0. 


For a* = +A > 0, the demographic noise £, leads to 
extinction in a finite time that diverges as cry -A 0 ( 3 ] [ 6 ]. 
The external noise p reduces the most probable value of 
the population size, that becomes very close to zero when 
a v > yj2A/b [51] , 

The population stability can be quantified by the mean 
time to extinction T averaged over realizations starting 
at u( 0). For Eq. (2.1) with D = 0, T is given by fITj . 


T = 


exp (/" T(it)du) 

W) 


dvdz. 


(4.1) 


where 'F(u) = 2 M(u)/V(u), with M{u) = au — b u 2 + 


a 2 u /2 and V(u) = a 2 u 2 + a 2 u. The results of Eq. (4.1) 


are in good accord with those from numerical simula¬ 
tions, as illustrated in Fig. [ 6 ] When the noise intensity 
decreases, the time to extinction always increases, being 
divergent in the limit try -A 0. 


A. From local to global behavior 

In this section we investigate the effects introduced by 
patch coupling, i.e., when 5 / 0. Nonlocal contributions 
redistribute the individuals in space, driven by density 
and quality gradients. In Fig. [7]we show that D / 0 pre¬ 
vents the extinction events that occur when D = 0 (see 
Fig# Therefore, in contrast to the deterministic case, 
now spatial coupling is constructive. On the other hand, 
noise has also a constructive role when d/ 0 , differently 
to the uncoupled case, not only preventing extinction but 
also contributing to the increase of the population (as in 
the case D = 10). In a previous work i7], we already ob¬ 
served the constructive role in population growth of lin¬ 
early multiplicative Stratonovich noise in contrast with 
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the destructive behavior of its Ito version. Therefore, en¬ 
vironmental noise and coupling have a positive feedback 
effect on population growth, as shown in Fig. [7] 

We will compute the long-time total population den¬ 
sity Moo = lim^oo 77(f), which is useful to be compared 
with the initial value ?7o = 77 ( 0 ) = pL 2 uo = pL 2 (A/b), 
that represents the asymptotic total density in the deter¬ 
ministic uncoupled case. Then we will measure the ratio 
E = (?7oo}/?7'0, that represents a kind of efficiency, where 
the brackets indicate average over landscapes and noise 
realizations. 


In the upper panel of Fig. [8] we plot the ratio E as 
a function of D. We see that for very small values of 
D, the population is non null, although the ratio E is 
smaller than one. Moreover, for given D , the ratio E is 
smaller when the diffusive component is absent (<5 = 0). 
In all cases the ratio first increases with D and even ex¬ 
ceeds the value E = 1, indicating again that not only the 
noise has a constructive role in preventing extinction but 
also in promoting the increase of the initial total popu¬ 
lation. When diffusion is present (<5 > 0), the increase of 
E occurs up to an optimal value of the coupling D (with 
E > 1), above which the ratio decays. Hence, there is a 
nonlinear effect that does not reflect the linear combina¬ 
tion in Eq. (2.2), as shown in the lower panel of Fig. [8] 
The diffusive component, despite being much less effi¬ 
cient, like in placing individuals in unfavorable regions, 
acts with greater connectivity. Then, for small D , the 
nonlocal contribution of the diffusive coupling is much 
higher than in the <5 = 0 case, leading to a higher popu¬ 
lation size. In fact, the abrupt transition in the connec¬ 
tivity of the spatial coupling is mirrored in the abrupt 
change suffered by the ratio E as 6 becomes non null. 
Contrarily, for high values of D, 6 = 0 is more efficient 
due to high damage caused by an intense dispersal to¬ 



Figure 7. Temporal evolution of n/no for 5 = 0.5, p = 0.1, 
4 = 0.5, p — 2.0, a v = <7£ = 1 and values of D indicated on 
the figure. We highlight a single realization (black full line) 
for each set of 50 realizations (gray lines). The dashed line at 
77 = no is plotted for comparison. 


wards unfavorable regions, which in the case of Fig. [8]are 
the majority of the sites. All these observations highlight 
the importance of the diffusive strategy, that can become 
more efficient than the ecological pressure driven by the 
quality gradient. 


B. Habitat topology and coupling range 


The nonlocal contribution results from the combina¬ 
tion of the spread strategies, interaction range and topol¬ 
ogy, characterized by <5, 4 and p, respectively. Fig. [9] 
shows the ratio E as function of p with different values 
of 4 for <5 = 1 and <5 = 0. 

E > 1 means that the combination of habitat topol¬ 
ogy and spatial coupling range leads the population to 
profit from the environment fluctuations, increasing its 
size. The region E > 1 is bigger when individuals are 
selective with respect to their destinations (5 = 0) and 
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Figure 8. Ratio E = {n^/no as a function of D (upper 
panel) for different values of <5, and A as a function of 5, for 
different values of D (lower panel), and p = 0.1, 4 = 0.5, 
p = 2.0 and a v = = 1. The symbols represent the average 

over 20 samples and the vertical bars the standard error. The 
dashed line at E = 1 is plotted for comparison. 
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Figure 9. Ratio E = (n 00 )/no, as a function of p for different 
values of l c , when 5 = 0 (upper panel) and 5 = 1 (lower 
panel), with p = 0.1, D = 20 and 07 = oj = 1. The symbols 
represent the average over 20 samples and the vertical bars 
the standard error. The dashed line at E = 1 is plotted for 
comparison. 


increases with t c . For the diffusive strategy (5 = 1), 
E > 1 is attained only in a clustered habitat (large p) 
together with short-range dispersal (small £ c ). We have 
already seen that in a sparse habitat, diffusion represents 
a waste, specially if the dispersal is long-range. Instead, 
when 5 = 0, the habitat does not need to be so clustered 
or the range so short for population growth. In this in¬ 
stance, the optimal combination occurs in a clustered 
habitat but with long-range coupling. Finally note that, 
as l c increases, E becomes independent of the topology. 


p not only changes the proportion of favorable patches 
but also reshapes the distribution of distances between 
favorable patches. In Fig. [To] we show three different 
outcomes of the spatial structure and the corresponding 
distance distribution for a fixed value of p = 2. For low /?, 
patches organize in a kind of archipelago structure, that 
is much smaller than the system size, and the distance 
resembles that obtained for large p when p = 0.1. For 
high p , many points of the domain are visited creating a 
distance distribution that approaches the homogeneous 
form. For p higher than the value of the figure, pro¬ 
files very similar to those shown in Fig. [lO] are obtained. 
Meanwhile, for small values of p , the distribution is al¬ 
most invariant with p , being very close to that of the 
uniform case. This is due to frequent flights with lengths 
of the order of system size. Concerning the factor 7 ^ 
that reflects the topology, as defined in Eq. (A3), it can 
be affected by p more through the amount of favorable 
patches n v than by its indirect consequences on the spa¬ 
tial distribution P^. 


p — 0.01 p — 0.1 p — 0.8 
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Figure 10. Spatial structure and probability distribution of 
the minimal distance between favorable patches (averaged 
over 100 landscapes), for p = 2 and three different values 
of p indicated on the figure. 


C. Density of favorable patches 

Another important issue is the influence of the density 
p of favorable patches in the dynamics. Until now, we 
have kept it constant to highlight the effects of the hetero¬ 
geneity of the habitat and of the coupling schemes in the 
longtime behavior of the total population size. In terms 
of the protocol used to generate the ecological landscape, 


In Fig. |TTJ we show the ratio E as a function of p for 
the case p = 2. By comparing the outcomes for different 
values of 5, we see the impact of distinct connectivities. 
In order to interpret this figure, recall that the initial 
population density no is proportional to the number of 
favorable patches n v , namely uq = n v A/b = pL 2 . 

For 5 = 1, E presents a minimum value for p — 0.15. 
Beyond this value, E grows with p attaining the value 

























of the full favorable lattice. In the opposite limit of van¬ 
ishing p (no favorable patches), E diverges as far as, ac¬ 
cording to the model, (intrinsically) favorable patches are 
not necessary to promote growth due to the noisy growth 
rate. However, if noise is reduced, then the stochastic 
dynamics approaches the deterministic one, where the 
population will certainly go extincted. 

Now, turning our attention to the (5 = 0 case, E is 
monotonically increasing with p 1 also attaining a limiting 
value when p —» 1. Differently from the diffusive case, 
there exists a critical value p c = 4 x 10 -4 (n v = 4) for 
population survival. 

For small p , it is curious that the role played by the 
connectivity, according to the model, makes the diffusive 
behavior more efficient, while selective moves are impor¬ 
tant at high values of p. In this case, when the sys¬ 
tem is approaching a fully favorable landscape, the ratio 
E tends to be the same for different values of 5. For 
intermediate values of p , we see that the selective strat¬ 
egy overcomes the diffusive one (but never overcomes the 
combined scheme). 



Figure 11. Ratio E = (n 00 )/no as a function of the favorable- 
patch density p, for L = 100. Different values of 5 were 
also considered as indicated on the figure. The remaining 
parameters are A = b = 1, D = 20, t c = 0.5, p = 2.0 and 
a v = <7£ = 1. The symbols represent the average over 20 sam¬ 
ples and the vertical bars the standard error. The horizontal 
line represents E = 1. 


V. CONCLUDING REMARKS 

We implemented a general model in which the local 
dynamics, ruled by the canonical model m , was cou¬ 
pled through different schemes on top of a complex land¬ 
scape. This setting allowed to study the role of the habi¬ 
tat spatial structure and the stochastic fluctuations on 
the long-time state of the metapopulation. We restricted 
the analysis to a region of parameter space relevant to 
display the main features and the interplay between the 
different processes involved. For the deterministic case, 


we have shown that, for small spread rates D , the dis¬ 
tribution of favorable patches must be clustered enough 
for survival, while below a critical value of p extinction 
occurs. For the stochastic case, we have shown that noise 
in combination with spatial coupling has a constructive 
role, that drives the population to survival, in contrast 
to the decoupled case where isolated patches would be 
extincted in finite time. We also studied the effects of 
the spreading strategy, pointing out that a mixed strat¬ 
egy (diffusive dispersion plus selective routes) will result 
in a larger population size (Fig. [8]). 

Furthermore, we analyzed the ratio E as a function 
of the coupling range £ c and landscape parameter p for 
different dispersal strategies. The more clustered, the 
more viable the environment is. For the selective strat¬ 
egy, the coupling range improves E , being more effective 
in disperse environments. The effect of the range satu¬ 
rates probably due to the rapid exponential decay of the 
weight function. In contrast, for the diffusive strategy, 
the coupling range plays and opposite role, as far as it 
drives individuals to unfavorable regions. 

Our model could be improved in several directions. For 
instance by considering correlated environment fluctua¬ 
tions, exhaustible resources, etc. But, despite simple, the 
model shows the impact of spatial coupling, spatiotem- 
poral fluctuations and their interplay, allowing to foresee 
the conditions for population survival as well as the op¬ 
timal dispersal strategy. 


Appendix A: Stability of deterministic steady states 

In order to study how steady state stability is affected 
by spatial coupling, let us assume that the population is 
located at the favorable patches, which is true for small 
D (that is, close to the uncoupled case), and that the cou¬ 
pling is purely diffusive (5 = 1). For a favorable patch, 
the deterministic form of Eq. ( |2.1[ ) reads 

Ui = Aiii -bui + D - Ui)y(dij) 

= (A- D)iii- bUi + D'^u j y(d i j) , (Al) 

recalling that l{dij) = 1. To estimate the last term, 
that represents the flow of individuals from the neighbor¬ 
hood towards patch i, J?", we consider that Uj ~ u,. In 
this case 

£ -.{//,,). (A2) 

where the sum effectively runs over the n v favorable 
patches. The average over arrangements of a landscape 
7 m = can be estimated as 

7 li = n v jPp(e)e~V t ‘de. (A3) 
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It depends on p and on the density p, such that it varies 
from p (when p = 0 ) to 1 , in the extreme cases of either 
maximal density or very large p. That is, 7 ^ in creases 
with p, with p and with t c too. Then, Eq. (All can be 
approximated by 


iq ~ (A - D[ 1 - 7 |U ])uj - bul = Gui - bv% . (A4) 


indicating decrease of the population. In fact notice in 
Fig. [5] that the smaller D the less frequent the extinction 
events for a given p. This effect can be mitigated by 
the landscape, through parameter 7 ^, when the density 
of favorable sites or clusterization associated with large 
p increases. Eq. |A5| also provides the linear stability 
condition for the null state. If G < 0, the population 
will decrease and go extincted. 


If G > 0, the population will grow and assume a finite 
value, bounded by the carrying capacity. Meanwhile, D 
diminishes the effective growth rate G, that becomes neg¬ 
ative for sufficiently large D , namely for 

D > A/( 1 - 7 ^) (A5) 
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